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ABSTRACT 

A fast ccwvuter program, GRID3C, has been developed to generate multilevel 
three-dimensional, C-type, periodic, boundary conforming grids for the calcu- 
lation of realistic turbomachinery and propeller flow fields. The technique 
is based on two analytic functions that conformally map a cascade of semi- 
infinite slits to a cascade of doubly Infinite strips on different Riemann 
sheets. Up to four consecutively refined three-oimensional grids can be auto- 
matically generated and permanently stored on four different computer tapes. 
Grid nonorthogonality is introduced by a separate coordinate shearing and 
stretching performed in each of three coordinate directions. The grids can be 
easily clustered closer to the blade surface, the trailing and leading edges 
and the nub or shroud regions by changing appropriate input parameters. Hub 
and duct (or outer free boundary) can have different ax i symmetric shapes. A 
vortex sheet of arbitrary thickness emanating smoothly from the blade trailing 
edge is generated automatically by GRID3C. Blade cross-sectional shape, chord 
length, twist angle, sweep angle, and dihedral angle can vary in an arbitrary 
smooth fashion in the spanwise direction. Input coordinates must be Cartesian, 
while the output grid coordinates can be Cartesian or cylindrical. 


INTRODUCTION 

When numerically solving partial differential equations governing the flow 
of fluid through realistically shaped configurations, exact boundary condi- 
tions must be applied at correct locations. This is especially important when 
calculating internal flows and flows that are governed by nonlinear partial 
differential equations. Seemingly negligible alterations of geometrical shape 
or flow conditions at the boundary can drastically change the basic features 
of the flow field, for example, choking an originally unchoked flow or chang- 
ing a shock-free flow into a shocked flow^. The most economical and accu- 
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rate way to numerically apply exact boundary conditions on solid boundaries is 
to generate and use a coo^utational grid that conforms to these surfaces (f'ig> 
1). Recent numerical techniques do not require orthogonal grids because 
they use locally isopararotric formulation when numerically determining deri- 
vatives of gecxnetric and flow variables. A widely accepted procedure for ac- 
celerating an iterative solution process of the flow equations and for resol- 
ving or capturing high flow gradients is to perform calculations on a sequence 
of several successively refined grids. The multigrid technique usually 
requires four to six such grids. For realistic threedinwnsional configura- 
tions the number of grid prints to be generated is prohibitively large even 
for inviscid flow calculations. Cofi^utational grids for such configurations 
should be easy to regenerate if shock waves and vortex sheets are to be better 
resolved or if the configuration of the solid boundaries changes with time. 

An H-type grid (fig. 1) provides excellent resolution of the flow field at 
upstream and downstream infinity. It is also the simplest grid to generate. 

At the same t‘me, H-type grid does not provide for an accurate treatment of 
rounded leading and trailing edges and wastes points in the flow domains away 
from the boundaries. An 0-type grid represents the other extreme. It gives a 
very poor resolution at infinities^, thus creating a problem when Cauchy- 
type boundary conditions must be enforced at the supersonic inflow boundary 
(fig. 2). A grid of the 0-type also does not provide desirable resolution in 
the vicinity of the vortex sheet. An open trailing edge simulation of the 
boundary layer displacement thickness effect cannot be readily incorporated. 
Nevertheless, an 0-type grid provides for accurate discretization of arbitrar- 
ily blunt leading and trailing edges and requires a minimum number of grid 
points. A combination of an 0-type grid in the upstream region and an H-type 
grid in the downstream region creates a C-type grid. This type of grid pro- 
vides for a good treatment of all boundary and periodicity conditions in- 
cluding wake treatment and supersonic exit flow, although it lacks an adequate 
resolution at upstream infinity (fig. 1). 

In turbomachinery and rotorcraft flow field calculations the flow field is 
periodic and a geometrically periodic grid provides for a simple and accurate 
way to enforce the flow periodicity. The simplest and fastest way to generate 
nonorthogonal periodic grids is to avoid time-consuming techniques based on 
the numerical solution of sets of partial differential equations whenever pos- 
sible. Instead, a basic knowledge of complex variables and conformal mapping 
can be used together with a few additional nonorthogonal coordinate shearings 
and stretchings. A three-dimensional, periodic, 0-type grid generator code 
was already developed^ by using this technique, which guarantees that the 
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grid lines of the same family do not intersect because the basis of the tech- 
nique is conformal mapping. Another view of a three-d1n«ns1onal, periodic 
0-type grid is presented in figure 3. 

The Computational Fluid Mechanics Branch of the NASA Lewis Research Center 
provided computational Facilities used in this work. Dr. Charles Putt of NASA 
Lewis Computer Services Division, Dr. Bharat Soni of Sverdrup Technology, Inc., 
and Mr. William Usab of MIT and United Technologies Research Center exercised 
the computer codes and provided several grid displays. 

SHEARING AND STRETCHING IN PHYSICAL SPACE 

Conformal mapping can be applied only to two-dimensional plane surface 
problems. A general procedure for creating such planes can be best described 
in the case of a rotor mounted on a hub shaped like a doubly infinite circular 
cylinder and confined inside a doubly infinite circular-cylinder-shaped duct. 

The intermediate doubly infinite circular-cylinder-shaped surfaces intersec- 
ting the blades can be viewed as planes when expressed in terms of (x,r«) co- 
ordinates. A standard procedure for creating three-dimensional blade shapes 
is to specify local airfoil shapes on a number of input planes that are or- 
thogonal to a straight radial line. This radial line (z axis in fig. 4) is 
called a stacking axis, and local blade sweep and dihedral angles are measured 
from that line (fig. 1). Input planes are defined by z = constant. Inter- 
mediate cylindrical surfaces, which we seek for the next step in this grid 
generation procedure are aefined by r = constant. To obtain an intersection 
contour between the blade surface and r ® constant cylindrical surfaces, a 
spline fitting and interpolation procedure is used along the blade. Input 
airfoil (x^,y^) coordinates on z = constant planes are transformed into 
cylindrical coordinates 



(1) 

= arc tan(y^ /z. ) 

(2) 

- (y^ 

(3) 


Cylindrical coordinates (x,re) are interpolated at r « constant spanwise loca- 
tions, thus defining blade cross sections on r * constant cylindrical surfaces. 

On the other hand, realistically shaped hubs and ducts are not doubly in- 
finite circular cylinders but axisymmetric surfaces. Therefore, the inter- 
mediate surfaces are also axisymmetric and not cylindrical. Nevertheless, the 
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same grid generation technique can be used if a sinple nonorthogonal shearing 
(or normalization) and stretching of the radial coordinate (fig. 3) Is per- 
formed. Nonorthogonal (unidirectional) shearing of the r coordinate converts 
the axisymraetric surfaces into cylindrical surfaces defined by R ■ constant. 
Let subscripts H,T, and D designate R » constant surfaces corresponding to 
hub, blade tip, and duct (or outer free boundary) location, respectively. Also 
let the normalized radial coordinate be defined as 

r (xj - ru(xJ 

The radial coordinate in the hub-to~t1p region is stretched and sheared with 
the following function 

R = + (R^ - R^)((R/R^) + A sin(2ir R/R^)) (5) 

The following value was obtained from experience 

R^ » N/50.0 (6) 

The stretching parameter. A, gives best results if 

0.12 > A > 0.0 (7) 

When A = 0, the cylindrical cutting surfaces R = constant are equidistantly 
spaced from hub to tip. Let the normalized, sheared radial coordinate in the 
region between the blade tip and the duct (or outer radial boundary) surface 
be 

R = (R - R^)/(Rq- R^) (8) 

The stretching function for the tip-to~duct domain is chosen to be 

R = 1.0 + (R^ - o)R + q R^ (9) 

This function must have the same slope, q, at the location R » 1 as the 
stretching function in the domain between the hub anc the tip (eq. 5). 

q * (1 + A)(l - R^)/R^ (10) 

Combining the two stretching functions (eqs. 5 and 9) provides for a smooth 
and continuous transformation from the physical r coordinate into the sheared 
R coordinate (fig. 4). For a stator or rotor with no tip clearance, equation 
9 is not needed. 
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Frequently, the input points are not clustered in the same regions on each 
input plane. Moreover, the number of input points defining the blade cross 
section on each input plane can vary from one input plane to the next. To 
accurately determine intersection contours between the blade surface and the 
axisymmetric surfaces, the corresponding input points must be located at the 
same percentage of the blade chord length on each input plane. Implicitly, 
this means that the number of input points must be the same on all input 
planes. Therefore, these input points must be appropriately redistributed on 
each input plane. This redistribution can be performed with respect to the 
input airfoil contour coordinate defined as 

S - [(»i - <11) 

Then the input Cartesian coordinates can be expressed in terms of the input 
airfoil contour coordinates. Coordinate s is measured clockwise around the 
input airfoil contour, starting and ending at the trailing edge point. As it 
was stated earlier, the number of contour points on the pressure surface must 
be the same as the number of contour points on the suction surface. For non- 
symmetric airfoils the lengths of these two contour lines are generally not 
the same. Let ITS denote the trailing edge point on the suction side and ITP 
denote the trailing eoge point on the pressure side of the input airfoil. 

Also, let LE denote the leading edge, that is, the point that is farthest from 
the trailing edge. The normalized surface coordinate is defined as 
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^LE ■ ^ITP 


(12) 


The redestribution of input points is performed with the following stretching 
function 


S = (1 - 5)5® + 5[1 - (1 - 5)®] (13) 

where the exponent B should satisfy 

1.4 > B > 1.0 (14) 


When B = 1 the points are equidistantly spaced along the airfoil contour. The 
points along the pressure surface are redistributed by using the formula 


S 



^ITP^ 


(15) 
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and the points along the suction surface are redistributed by using the formula 
^ ~ ^IIP^ 

This redistribution of input coordinates x and y is performed with a cubic 
spline fitting applied in the s direction. Interpolation is performed at S 
locations. Spline fitting and interpolation are also used with respect to the 
R coordinate in order to find the points on intersection contours between the 
blade surface and the intermediate axisymmetric surfaces. Locations of those 
points in the physical space will not be altered with the subsequent mapping- 
remapping procedure. 

The exact shape of the wake of arbitrary thickness is not known a priori. 

To eliminate the need for specifying the location of the wake in the prepara- 
tion of the input, the shape of the wake centerline is automatically generated 
by using the simple polynomial expression 


y = a(x b{x - x^^)^ + c(x - x^^) + y^^ (17) 

Here the trailing edge point coordinates are 

’'tE = ("ITP * 
and 

The point where the wake centerline intersects the downstream-infinity cutoff 
boundary is defined witli the subscript EX. Let c be the average slope of the 
pressure and suction surfaces of the airfoil at the trailing edge , and let d 
be the slope of the expected flow angle at the exit boundary. Then the con- 
stants a and b in equation 17 are 


a = [x (c d) - 2 y ]/x^ 
and 

b = [3y - X (2c + d)]/x^ 

where 

\ = ^EX ■ *TE 


( 20 ) 


( 21 ) 


( 22 ) 
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and 

“ ^£X ■ ^TE 

Wake surface grid points are redistributed (stretched) with the formula 
X* « (x - Xjg)/x^ - n sin(ir(x - x-j-^)/x^) (24) 

The stretching exponent, n, is determined from the continuity of the slope of 
the stretching functions at the trailing edge (eqs. 13, 15, 16, 22, and 24) 

n = 1.05(1.0 - B/2x^)/it (25) 

w 

If the wake has a finite thickness, that is, if the blade trailing edge is 
open, coordinates of the points on the upper and lower surfaces of the wake 
are determined by adding and subtracting the trailing edge half thickness. 

The axial coordinate of the upper surface of the wake is determined from the 
formula 

x'^ = X + (Xjj^- Xjjp)/2 (26) 

and that of the lower surface of the wake by the formula 

X = X - (xjj^— Xjjp)/2 (27) 

with similar expressions for the y coordinate. Superscripts u and 1 designate 
the upper and lower surfaces of the wake, respectively. 

CONFORMAL MAPPING AND REMAPPING 

The conformal mapping portion of the present procedure for generating 
three-dimensional, periodic C-type grids was originally used by Sockol to 
generate orthogonal, two-dimensional, cascade C-type grids. If the blades are 
straight, semiinfinite twisted plates of zero thickness, their intersections 
with circular cylinders generates doubly infinite cascades of semi-infinite 
straight slits on each of the (x,Rg) planes (fig. 5). Each of these R « con- 
stant planes can be defined in terms of complex variables 

w * X iRe (28) 

The goal is to generate a boundary-conforming, periodic C-type grid on each of 
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the planes. Tni ■ task is iccomplished by conformally mapping the w plane via 
an intermediate "circle" duplex plane (fig. 6) 

V = c in (29) 

into the interior of a "nouuly infinite strip" plane (fig. 7) 

u = X + iY (30) 

Uniform grid in tne u plane is then conformally remapped into the w plane, 
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thus generating ihe desired C-type grid. As shown by Sockol a single ana- 
lytic function 

w = + e ^ -|l-(2s sin 6 2 cos B ln(2 cos $)) 

" e’‘®(i:i V - in ) - 2 cos e In (1 - v) (31) 

whet'.^ N is the number of blades and e is the local stagger angle on the R « 

constant surface, conformally maps the interior of the unit circle in the v 
plane to the interior of a periodic strip enveloping a semi-infinite slit in 
the w plane. The center of the circle (v = 0) maps into upstream infinity in 
the w plane and the point v = -1 maps into downstream infinity in the w plane. 
The zero-thickness slit between the points v = 0 and v = -1 maps into the up- 
per and lower periodic boundary of a periodic strip in the w plane. The 
circle in the v plane maps into a semi-infinite straight slit in the w plane. 
A doubly infinite cascade of semi-infinite straight slits in the w plane is 
thus created by conformally mapping a doubly infinite cascade of Riemann 
sheets (v planes) that are interconnected through the slits between the points 
V = 0 end V = -1. Sockol^ used a simple analytic function 

V = tanh(u^/2) (32) 

to conformally map the interior of a doubly infinite straight strip in the u 
plane into the interior of a unit circle in the v plane. The lower strip 
boundary (Y = -iir /2) in the u plane maps into the circle in the v plane. 
The upper strip boundary (Y = 0) maps into a zero-thickness slit between the 
points V = 0 and v = -1. Axial infinities (X = ±®) map into a single point 
(v = -1). The origin (X = 0;Y = 0) in the u plane maps into the origin (v » 
0) in the v plane. 
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Realistically shaped blade airfoils are not straight semi-infinite lines 
of zero thickness. A C-type grid generated with the use of equations 31 and 
32 alone will not conform to the actual airfoil cascade shapes on R - constant 
surfaces. To generate a C-type grid that conforms to the shape of the airfoil 
and wake , several nonorthogonal coordinate shearings and stretchings are used. 
Airfoil surface points are conformally mapped from the w plane via the v 
plane into the u plane. As a result, the circle in the v plane becomes de- 
formed (fig. 7), and the corresponding lower wall in the u plane becomes an 
irregular line (fig. 8). The inverse of equation 31 cannot be analytically 
obtained for staggered cascades. Therefore, a Newton-Raphson procedure is 
used to iteratively evaluate on a point-by-point basis the pairs of (t»n) 
coordinates corresponding to the given pairs of (x,R«) coordinates. By using 
an analytic inverse of equation 32, that is, 

U . [1„ (33) 

the deformed circle is conformally mapped from the v plane into the u 
plane. 

SHEARING AND STRETCHING IN COMPUTATIONAL SPACE 

It should be pointed out that with the increase in stagger angle in the w 
plane the image of the leading edge point shifts along the deformed circle in 
the V plane and along the deformed lower boundary in the u plane. To in- 
sure that the corresponding points along the periodic boundaries in the w 
plane have the common values of x coordinate, their images in the u plane are 
placed symmetrically along the Y * 0 line (fig. 9). At the same time these 
periodic points are distributed with a simple stretching function 

- e sin(2n X^(X^.^^ - x^^p)) (34) 

Superscript U denotes the upper wall (Y = 0) of the u plane and superscript 
L denotes the lower irregular boundary of the u plane. The stretching 
coefficient e is decermined from experience as 

e = 0.18 - 0.05 ln(2Rii/Nt) (35) 

where t is the local blade chord. The periodic grid points located in the 
wake region are redistributed by using the expression 


9 


x'J . x'^ - f sin(2.('Xl‘^- Xnsl'I'^MAXXP’ *ITS>> 

where MAXXP denotes the last point on the upper surface of the wake. The 
stretching coefficient, f, is detennined also from the experience as 

0.10 > f > 0.0b (37) 

Because only a finite length of the wake is conformally mapped from the w 
plane into the u plane, the deformed strip in the u plane has a finite 
length. The shape of the end wall boundaries in the u plane are determined so 
that they meet the lower boundary of the strip in the u plane almost orthogo- 
nally (fig. 8). Consequently, grid orthogonality is well preserved at the 
wake. Coordinates of the grid points inside the strip in the u plane are de- 
termined from 


Y = y‘-((Y/y'-) + g sin(itY/Y'-)) (38) 

and 

X = X^ + (X*- - X^) ((Y/Y^) + C sin(ttY/Y'-)) (39) 

where 

0.30 > C > 0.16 (40) 

g = C (1.0 -1.0/cosh h) (41) 

h = b (X>J/X„%p ) (x,2) 


Stretching coefficients C, g, and h are determined from experience and from 
the condition that C-type grid lines in the w plane closely follow the wake 
contour. Larger values or C generate grids suitable for viscous flow calcula- 
tions, because grid layers are positioned closer to the blade and wake 
surface. 

The resulting two-dimensional nonorthogonal periodic grid in the u plane 
is conformally mapped back into the w plane on a point-by-point basis. Final- 
ly, determination of the physical r coordinates of the grid points on the 
(x,Re) planes is obtained by reshearing the n coordinate (eqs. 4, 5, 8, and 9) 
and fitting it with respect to the x coordinate with a cubic spline. 


RESULTS 

On the basis sf the preceding analysis, a con^uter program GRID3C was de- 
veloped and tested^. Program 6R103C consists of USO card statements and 
requires approximately SOO K of computer memory Because of the analytical 
character of most of the transformations used, GR1D3C is very fast. To gener- 
ate and permanently store x,y,z coordinates of a typical four-level grid se- 
quence consisting of (33*8*6), (63*13*11), (123*23*21), (243*43*41) grid points, 
respectively, GR103C requires between three and four minutes of CPU time on an 
IBM 370/3033 coiT4)uter. The Newton- Raphson iterative point-by-point mapping 
procedure of the airfoil and wake contour from the w plane Into the v plane 
consumes most of the computer time. But this procedure needs to be performed 
only once on each axisymmetric surface. 

Input to GRID3C must be provided in the x,y,z coordinate system, while the 
output grid coordinates can be computed In the x,y,z or x,r,e coordinate 
system. GRI03C can automatically generate up to four successively refined 
three-dimensional grids and store them on four separate tapes. Computational 
grids for che blades with closed trailing edge (fig. 10) and for the blades 
with open trailing edge (fig. 11) can be generated with 6R1D3C code. For re- 
petitive runs with different numbers of blades or different blade setting an- 
gles, only one input parameter needs to be changed in the Input deck. Clus- 
tering of grid points closer to the leading and trailing edges and closer to 
the blade and vortex sheet surface (fig. 12) can be easily achieved by varying 
coordinate stretching parameters A, B, and C. Grid nonorthogonality Is almost 
entirely removed from the airfoil and wake surface. Nevertheless, grid nonor- 
thogonality can become intolerable if this grid generation technique is ap- 
plied to closely spaced, highly staggered and cambered blades. Nonorthogo- 
nality can become excessive in the leading edge region of any blade if the end 
point of the semi-infinite slit in '.he w plane is not positioned approximately 

midway between tne leading edge and Its center of curvature. 

An unsatisf actot y grid resolution inherent to the C-type grids can be ob- 
served in figure 13. This figure siiows a rectangular wing - cylindrical fuse- 

lage combination and two computational grid surfaces: one corresponding to the 
surface of the fuselage and the other being in intermediate surface located 
between the hub and the wing tip. Note that the wing extension beyond the tip 
has linearly increasing cord length. The GRI3C code automatically calculates 
wing (or blade) chord lengths at the off-tip locations based on the constraint 
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that gap-to-chord ratio at the tip should be retained at all outer spanwise 
locations. Key elements of a three-dimensional C-type grid generated by the 
GR1D3C code for an advanced, eight-blade, transonic, NASA propeller Is 
presented In figures 14 and 15 with Intersection contours between a blade and 
the axl symmetric sur- faces shown. Note the large twist, sit^ep, and taper 
variations and the fact that the propeller hub Is axisymmetric. 

With r.tinor modifications GRI03C can be used for generating con^utational 
grids applicable to a midmounted wing-body combination or a finned missile In 
free air or Inside a wind tunnel having axisymmetric walls. 
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Figure 1. - Three basic types of two-dimensional, conforming, computational 
grids. 
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Figure 2 . - Axisynmetric view of a three-dimensional, 0-type, periodic bound 
ary conforming grid for NASA eight-blade transonic prop fan. Shown are the 
hub surface grid and three neighboring blades with their surface grids. 
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Figure 7 
Deformed 
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Figure 8. - "Strip" plane obtained by conformally mapping "circle" plane. 
Upper boundary corresponds to periodic boundaries, and lower boundary to air- 
foil shape. 
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Figure 9. - Nonorthogonal coordinate shearing and stretching concept applied 
to X (eqs. 34, 36, and 39) and Y (eq. 38) coordinates results in a desired 
rectangular computational surface. 
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Figure 10. - f^t\ example of a two-dimensional (x,Re) surface discretized with a 
coarse, C-type, periodic grid. 


23 


OR»aiNM-MjW* 

OFPOORQUW^ 



Figure 11. - Two-dimensional (x.Re), C-type, periodic, 
grid for a cascade of blades with open trailing edges. 


boundary conforming 
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Figure 12. - Effect of controlled grid clustering. Grid points can be easily 
concentrated in the regions of leading and trailing edges as well as closer to 
the surface of the airfoil and its wake. 
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Figurf 13. - Elements of a three-dimensional, C-type, perioaic grig generated 
by GR103C code for a geometry consisting of a rectangular unswept wing attach- 
ed to a circular cylinder. Note deteriorating grid quality in the far Uj- 
stream region. Only every fourth cylindrical surface is shown. 
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Figure 14. - Blade surface grid and one of the axi symmetric surfaces generated 
by 6RI03C for an advanced, eight-blade NASA prop fan. 





